/*************************************************************/
/*    Source for the stress of a neo-hookean finite strain   */
/*    Initially based on ComputeFiniteStrainElasticStress    */
/*                  Status: not working                      */
/*                                                           */
/* Copyright (C) 2018  eDgar                                 */
/*                                                           */
/* This program is free software: you can redistribute it    */
/* and/or modify it under the terms of the GNU General       */
/* Public License as published by the Free Software          */
/* Foundation, either version 3 of the License, or (at your  */
/* option) any later version.                                */
/*                                                           */
/* This program is distributed in the hope that it will be   */
/* useful, but WITHOUT ANY WARRANTY; without even the        */
/* implied warranty of MERCHANTABILITY or FITNESS FOR A      */
/* PARTICULAR PURPOSE.  See the GNU General Public License   */
/* for more details.                                         */
/*                                                           */
/* You should have received a copy of the GNU General Public */
/* License along with this program.  If not, see             */
/* <http://www.gnu.org/licenses/>.                           */
/*************************************************************/

#include "ComputeStressNeohookean.h"

registerMooseObject("spiderApp", ComputeStressNeohookean);


// Just an idea to stop copying the code from ComputeElasticityTensorNeohookean
// template <>
// InputParameters
// validParams<ElementAverageValue>()
// {
//   InputParameters params = validParams<ElementIntegralVariablePostprocessor>();
//   return params;
// }


template <>
InputParameters
validParams<ComputeStressNeohookean>()
{
  InputParameters params = validParams<ComputeStressBase>();
  params.addClassDescription
    ("Compute neo-Hookean stress using elasticity for "
     "finite strains");

  ///////////////////////////////////////
  // Repeated code from                //
  // ComputeElasticityTensorNeohookean //
  ///////////////////////////////////////

  // Give a default value to \mu (set -1; if the value is
  // not set to a positive value by the user, there is a
  // check below which throws an error)
  params.addParam<Real>
    ("mu",
     -1.0,
     "Lamé second constant (shear) of the material.");
  // Allow lambda to be an input parameter
  params.addParam<Real>
    ("lambda",
     "Lame's first constant (non-shear) of the material.");

  // Allow Poisson's ratio and Young's modulus as input
  // parameters
  params.addParam<Real>
    ("nu", "Poisson's ratio of the material.");
  params.addParam<Real>
    ("E", "Young's modulus of the material.");

  return params;

  ///////////////////////////////////////
  // -- end Repeated code              //
  ///////////////////////////////////////
}

ComputeStressNeohookean::ComputeStressNeohookean(
    const InputParameters & parameters)
  : ComputeStressBase(parameters),
    GuaranteeConsumer(this),

    ///////////////////////////////////////
    // Repeated code from                //
    // ComputeElasticityTensorNeohookean //
    ///////////////////////////////////////

    // Get whether or not the parameters were given
    _lambda_set(parameters.isParamValid("lambda")),
    _mu_set(parameters.isParamValid("mu")),
    _nu_set(parameters.isParamValid("nu")),
    _E_set(parameters.isParamValid("E")),

    // Take the Lamé constants from
    // the GetPot (.i) file
    _lambda(_lambda_set ? getParam<Real>("lambda") : -1),
    // The same for `mu'
    _mu(_mu_set ? getParam<Real>("mu") : -1),

    // `nu' and `E' are parameters, because I don't need to
    // store them for anything else (so far).
    _nu(_nu_set ? getParam<Real>("nu") : -1),
    _E(_E_set ? getParam<Real>("E") : -1),

    // Take the deformation gradient from
    // ComputeFiniteStrain
    _deformation_gradient(getMaterialProperty<RankTwoTensor>
                          ("deformation_gradient")),

    ///////////////////////////////////////
    // -- end Repeated code              //
    ///////////////////////////////////////

    // "You can declare material property for J_1 in the
    // ComputeElasticityTensorNeohookean.C and get material
    // property J_1 in ComputeStressNeohookean.C." --Wen Jiang
    //
    // I hope that this is the way to use a property that I
    // declared in ComputeElasticityTensorNeohookean
    _J_1(getMaterialProperty<Real>
         ("jacobian_inverse")),
    _logJ(getMaterialProperty<Real>
          ("jacobian_log")),

    // See note about "Rotate the stress" in source code
    // (following ComputeFiniteStrainElasticStress.C and
    // ComputeFiniteStrain::computeQpIncrements).
    _rotation_increment(getMaterialProperty<RankTwoTensor>
                        ("rotation_increment"))
{
  ///////////////////////////////////////
  // Repeated code from                //
  // ComputeElasticityTensorNeohookean //
  ///////////////////////////////////////
  // Neo-hookean material model

  // Now we are going to initialise variables for the analysis
  // TODO:

  // Make sure that only one of these are given:
  // 1. Young's modulus and the Poisson's ratio
  // 2. Lamé's constants lambda and mu
  if (! ( (_E_set && _nu_set) || (_lambda_set && _mu_set) ) )
    mooseError
      ("One of these combinations is needed:\n"
       "(1) Young's modulus (E) and the Poisson's ratio (nu)\n"
       "(2) Lamé's constants lambda and mu\n"
       "in the input file for the neo-Hoookean\n"
       "'" + name() + "'.");
  else if (_E_set && _nu_set) {
      // If both Young's modulus and Poisson's ratio are
      // given, use them to calculate Lamé's constants
      if (_nu_set && (_nu <= -1.0 || _nu >= 0.5))
        mooseError
          ("Poissons ratio must be greater than -1 and less"
           " than 0.5 in material '" + name() + "'.");
      if (_E_set && _E <= 0.0)
        mooseError
          ("Youngs modulus must be positive in material '"
           + name() + "'.");

      // Value of Lamé constants in terms of Young's modulus
      // and Poisson's ratio
      _mu = _E / (2 * (_nu + 1));
      // // 3 multiplications, 1 division, 1 addition, 1 subtraction:
      // _lambda = E * _nu/ ((_nu + 1) * (1 - 2 * _nu));
      // 3 multiplications, 1 division, 1 subtraction, less clear:
      _lambda = 2 * _mu * _nu / (1 - 2 * _nu);
    }
  else if (_lambda_set && _mu_set) {
    // Make sure that the values of Lamé's constants are
    // not both negative or zero.
    // Holzapfel, G. A., Nonlinear solid mechanics : a
    // continuum approach for engineering (2000), :
    // Wiley. p.250
    if ((_lambda < 0.0 || _mu < 0.0) ||
        (_lambda == 0.0 && _mu == 0.0))
      mooseError("The Lamé constants must be positive in material '"
                 + name() + "'.");
  }
  ///////////////////////////////////////
  // -- end Repeated code              //
  ///////////////////////////////////////
}

void
ComputeStressNeohookean::initialSetup()
{
  if (!hasGuaranteedMaterialProperty
      (_elasticity_tensor_name, Guarantee::ISOTROPIC))
    mooseError("ComputeStressNeohookean can only be used "
               "with elasticity tensor materials "
               "that guarantee isotropic tensors.");
}

void
ComputeStressNeohookean::initQpStatefulProperties()
{
  ComputeStressBase::initQpStatefulProperties();
}

void
ComputeStressNeohookean::computeQpStress()
{
  // \LaTeX:
  /** \begin{align}
   *    \sigma_{ij} = \frac{1}{J}\left(
   *      \lambda_{0}\cdot\ln{\left(J\right)}\cdot\delta_{ij}
   *       + \mu_{0}\cdot \left(B_{ij}
   *          - \delta_{ij}\right)
   *    \right)
   *  \end{align}
   */
  // Identity matrix
  // 1  0  0
  // 0  1  0
  // 0  0  1
  // initIdentity() in
  // [[moose/framework/src/utils/RankTwoTensor.C]]
  const RankTwoTensor dij(RankTwoTensor::initIdentity);

  // // This works, but I don't want to recalculate, thus
  // // get it from ComputeElasticityTensorNeohookean.C)
  // // Jacobian (determinant of the deformation gradient)
  // Real _J = _deformation_gradient[_qp].det();
  // // 1 / _J
  // Real _J_1 = 1 / _J;
  //
  // Deformation gradient tensor F = 1 + du_{i}/dX_{j}
  const RankTwoTensor F = _deformation_gradient[_qp];
  // Inverse of the determinant of F (from ComputeElasticityTensor)
  const Real J_1 = _J_1[_qp];
  // Log of the determinant of F (from ComputeElasticityTensor)
  const Real logJ = _logJ[_qp];

  // Left Green-Cauchy deformation tensor
  const RankTwoTensor Bij = F * F.transpose();

  _stress[_qp] = J_1 * (_lambda * logJ * dij + _mu * (Bij - dij));

  // Rotate the stress state to the current configuration
  // (following ComputeFiniteStrainElasticStress.C and
  // ComputeFiniteStrain::computeQpIncrements).
  // We assume that the calculations of stress are done
  // without the rotation portion of F = R·U, and _stress
  // actually needs it. It would work awfully without it for
  // linear behaviour and not work at all for non-linear
  // materials.
  _stress[_qp] = _rotation_increment[_qp] \
    * _stress[_qp] * _rotation_increment[_qp].transpose();
  // _stress[_qp] = RankTwoTensor(RankTwoTensor::initIdentity) * 2;

  // "Jacobian_mult is a fourth order tensor that is the
  // partial derivative of your stress with respect to the
  // strain" --Michael Tonks
  //
  // The question was answered by Dr. Lujan with the help of
  // [[cite:BonetWood99][Eq. 5.13]]: \[\dot{\sigma_{ij}} =
  // c_{ijkl} : \dot{d_{kl}}\], where \(d_{kl}\) is the
  // push-forward of \(\dot{E}\), which implies that
  // \(c_{ijkl} = \frac{\dot{\sigma_{ij}}}{\dot{\d_{kl}}}\)
  //
  // Compute dstress_dstrain
  _Jacobian_mult[_qp] = _elasticity_tensor[_qp];
  // Also see: FiniteStrainHyperElasticViscoPlastic.C
  // _Jacobian_mult[_qp] = _tan_mod * _dfe_df * _df_dstretch_inc;

  // // WAIT: testing
  // // From ComputeFiniteStrainElasticStress.C, (and defined
  // // in ComputeFiniteStrain.C)
  // // Assign value for elastic strain, which is equal to the
  // // mechanical strain
  // _elastic_strain[_qp] = _mechanical_strain[_qp];
}
